home *** CD-ROM | disk | FTP | other *** search
/ Aminet 22 / Aminet 22 (1997)(GTI - Schatztruhe)[!][Dec 1997].iso / Aminet / mus / misc / MusicIn.lha / musicin / FastFT_float.c < prev    next >
C/C++ Source or Header  |  1997-09-23  |  2KB  |  95 lines
  1. /*------------------------------------------------------------------------------
  2.  
  3.     File    :    FastFT_float.c
  4.  
  5.     Author  :    Henryk Richter
  6.                 (based on Sources by Stéphane TAVENARD)
  7.  
  8.     $VER:   FFT_float.c  0.1  (17/09/1997)
  9.  
  10.     FFT support routines for double precision floating point values
  11.  
  12. ------------------------------------------------------------------------------*/
  13.  
  14. #define FFT_Buggs
  15.  
  16. #include <stdio.h>
  17. #include <math.h>
  18.  
  19. #ifdef FFT_Buggs
  20. #include "FastFT_float_asm.h"
  21. #endif
  22.  
  23. #include "FastFT_float.h"
  24.  
  25. static WORD INDEX[ FFT_RANGE_MAX ];
  26. static FLOAT SINCOS [ FFT_RANGE_MAX  ]; // FFT_RANGE_MAX / 2 * 2
  27. static FLOAT SINCOS2 [ FFT_RANGE_MAX  ];
  28.  
  29. static int range = 0;
  30. static int power = 0;
  31.  
  32. static void build_index( int power )
  33. {
  34.    WORD n, i, j, l, m, x;
  35.  
  36.    x = power - 1;
  37.    n = 1<<(x);
  38.    for( i=0; i<n; i++ ) {
  39.       INDEX[ i ] = 0;
  40.       l = 1;
  41.       m = n;
  42.       for( j=0; j<x; j++ ) {
  43.      m = m>>1;
  44.      if( i & l ) INDEX[ i ] += m;
  45.      l = l<<1;
  46.       }
  47.    }
  48. }
  49.  
  50. static void build_sin_cos( int power )
  51. {
  52.    WORD i, n, j;
  53.    double p;
  54.  
  55.    n = 1<<(power-1);
  56.    p = PI;
  57.    for( i=0; i<n; i++ ) {
  58.       j = INDEX [i];
  59.       SINCOS[ 2*i ]     = (float)sin( 2 * p * (double)j / (double)n);
  60.       SINCOS[ 2*i + 1 ] = (float)cos( 2 * p * (double)j / (double)n);
  61.  
  62.       SINCOS2[ 2*i ]     = (float)sin( p * (double)i / (double)n);
  63.       SINCOS2[ 2*i + 1 ] = (float)cos( p * (double)i / (double)n);
  64.    }
  65. }
  66.  
  67. void FastFT_FLOAT_forward( FLOAT *x_real, FLOAT *x_imag,
  68.             FLOAT *energy, FLOAT *phi, int n )
  69. {
  70.    if( n > FFT_RANGE_MAX ) {
  71.       fprintf( stderr, "FFT_DOUBLE_forward: n is out of range (%ld>%ld)\n",
  72.                n, FFT_RANGE_MAX );
  73.       return;
  74.    }
  75.    if( n != range ) {
  76.       range = 1;
  77.       power = 0;
  78.       while( range < n ) {
  79.      range = range << 1;
  80.      power++;
  81.       }
  82.       if( range != n ) {
  83.      range = 0;
  84.      power = 0;
  85.      fprintf( stderr, "FFT_DOUBLE_forward: n is not a power of 2 (%ld)\n", n );
  86.      return;
  87.       }
  88.       build_index( power );
  89.       build_sin_cos( power );
  90.    }
  91.  
  92.    ASM_FastFT_main_loop( x_real, x_imag, SINCOS, SINCOS2, power );
  93.    ASM_FastFT_energy_phi( x_real, x_imag, energy, phi, n );
  94. }
  95.